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INTRODUCTION 


The  purpose  of  this  research  was  to  develop  efficient  and  accurate  finite 
elements  for  the  large  displacement,  transient  analysis  of  shells.  It  was 
apparent  from  the  beginning  of  this  research  program  that  in  order  to  achieve 
these  goals,  minimal  quadrature  rules  would  have  to  be  developed  so  that  the 
fewest  possible  number  of  quadrature  points  would  be  used  in  an  element.  In 
the  case  of  the  4-node  shell  element,  minimal  quadrature  consists  a  single 
quadrature  point  per  element;  in  the  9-node  Lagrange  shell  element,  minimal 
quadrature  is  a  2x2  Gauss  quadrature.  Compared  to  full  quadrature,  reduced 
quadrature  reduces  solution  time  by  50%  to  75%  and  also  enhances  the  accuracy 
of  shell  elements. 

However,  minimal  quadrature  elements  have  one  important  shortcoming:  they 
possess  spurious  singular  modes,  often  known  as  hourglass  modes,  which  can 
completely  destroy  a  solution.  Therefore,  procedures  were  developed  for 
controlling  these  spurious  modes.  These  methods  have  been  called  y-methods 
and  they  involve  a  special  projection  so  that  the  consistency,  that  is  the 
ability  to  meet  the  patch  test,  of  the  finite  element  is  not  lost.  We  have 
been  able  to  apply  these  methods  to  both  linear  and  nonlinear  problems  as 
evidenced  by  the  results  reported  in  Belytschko,  Tsay  and  Lin  (1981).  Some  of 
the  nonlinear  results  obtained  in  that  paper  are  reviewed  in  Appendix  C  of 
this  report.  It  can  be  seen  from  these  results,  that  the  method  has  indeed 
become  effective  in  providing  accurate  nonlinear  solutions.  The  method  has 
already  been  Incorporated  in  the  general  purpose  program  ABAQUS. 

In  addition  to  the  basic  development  of  the  stabilization  procedure,  the 
major  findings  of  this  project  are:  (1)  the  identification  of  the  membrane 
locking  phenomenon  which  impedes  the  convergence  of  any  fully  integrated 
curved  element;  (2)  the  development  of  general  methods  for  ameliorating 


membrane  locking  through  both  explicit  mode  decomposition  projections  and 
through  implicit  projections  by  means  of  reduced  integration;  (3)  the 
development  of  stabilization  procedures  for  higher  order  elements  such  as  the 
9-node  element  which  satisfy  basic  consistency  and  the  patch  test. 

The  motivation  for  the  research  into  the  curved  elements  can  be 
understood  by  examining  Figure  la.  These  results  were  obtained  with  the 
widely-used,  9-node  Lagrange  C°-shell  element.  As  can  be  seen,  it  is 
difficult  to  choose  a  quadrature  scheme  for  the  9-node  element  which  is  both 
accurate  and  stable.  For  full  Integration  (3x3)  or  selective  reduced 
integration,  even  a  relatively  fine  mesh  such  as  this  results  in  errors  which 
are  unacceptably  large.  On  the  other  hand,  uniform  reduced  integration,  which 
in  this  case  is  2x2  on  all  terms,  provides  good  accuracy  but  results  in 
singularity  of  the  assembled  stiffness  for  some  support  conditions.  This 
limitation  of  uniform  reduced  Integration  probably  makes  it  unacceptable  for 
general  purpose  programs,  but  its  superior  accuracy  is  attractive. 

The  relatively  poor  performance  of  the  fully-integrated  (3x3)  9-node 
Lagrange  element  in  the  simple  arch  problem  is  actually  a  mild  case  of 
misbehavior.  As  shown  in  Fig.  lb,  for  more  complex,  deeply-curved  shells,  the 
behavior  of  both  the  3x3  and  selective-reduced  integration  schemes  can  be 
simply  abysmal .  In  this  case,  1445  degrees  of  freedom  for  a  quadrant  of  a 
shell  yielded  results  which  were  only  26%  of  the  exact  solution. 

A  similar  impasse  has  evolved  in  the  development  of  3-node,  18  degree-of- 
freedom  triangular  shell  elements.  Essentially,  prior  to  this  research  no 
element  of  this  genre  existed  which  could  solve  a  wide  class  of  shell  problems 
with  acceptable  accuracy.  While  certain  elements  performed  well  for  specific 
problems,  invariably  when  tested  on  a  set  of  demanding  shell  test  problems, 
their  performance  is  unsatisfactory. 
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Two  of  the  phenomena  which  have  been  identified  as  the  culprits  in  poor 
element  performance  are  shear  and  membrane  locking.  Shear  locking  was  first 
identified  by  Doherty  et  al .  (1969)  and  Ziekiewicz  et  al .  (1971),  who  at  the 
same  time  proposed  the  daring  remedy  of  reduced  integration.  Reduced 
integration,  and  its  offspring,  selective  reduced  integration  (SRI),  proved  to 
be  very  successful  in  ameliorating  the  effects  of  shear  locking  in  the 
analysis  of  plates  and  shells. 

The  term  "membrane  locking"  was  coined  by  Stolarski  and  Belytschko  (1981) 

(see  also  Stolarski  and  Belytschko  (1982)),  who  showed  that  it  is  related  to 
an  inadequate  representation  of  inextensional  deformation.  The  poor 
performance  of  many  elements  in  analyzing  the  response  of  shells  where 
inextensional  modes  of  deformation  are  important  has  been  noted  by  numerous 
authors,  including  Ashwell  and  Sabir  (1971),  Sabir  and  Ashwell  (1971),  Morley 
(1972),  Sabir  and  Lock  (1973),  Fried  (1973)  and  Dawe  (1974).  Problems  related 
to  inextensional  bending  were  also  discussed  later  by  Ashwell  (1976),  Lee  and 
Pian  (1978),  Noor  and  Peters  (1981),  Kikuchi  (1982),  Allman  (1982),  MacNeal 
(1982),  Morley  (1983)  Stolarski  and  Belytschko  (1983)  and  Kikuchi  and  Aizawa 
(1984).  Reduced  integration  is  also  effective  in  mitigating  membrane  locking; 
see  Parisch  (1979)  and  Stolarski  and  Belytschko  (1983). 

However,  prior  to  the  work  performed  in  this  contract,  membrane  locking 
in  complex  shell  elements  was  poorly  understood  and  its  pervasiveness  not 
appreciated.  Thus  many  users  of  degenerated  elements  did  not  anticipate  its 
appearance  in  those  elements,  for  the  idea  expressed  in  Noor  and  Peters 
(1981),  that  this  locking  "is  a  result  of  inadequate  representation  of  rigid 
body  modes"  was  widely  held.  Yet  degenerated  i soparametrics,  which  correctly 
incorporate  rigid  body  motion,  as  first  noted  by  Argyris  and  Scharpf  (1968), 
also  encounter  severe  membrane  locking. 


,>  f 


This  contract  has  shown  that  locking  can  be  eliminated  by  projection 
methods.  Among  these  are  mixed  formulations,  in  which  separate  interpolants 
are  used  for  the  stresses  and  displacements  as  in  Lee  and  Pi  an  (1978)  and  Noor 
and  Peters  (1981),  or  by  mode  decomposition  methods  in  which  the  nodal 
displacements  are  projected  so  as  to  minimize  parasitic  stresses.  Examples  of 
the  latter  are  the  Hughes  and  Tezduyar  (1981)  quadrilateral  plate  element  and 
the  Belytschko,  et  al .  (1984c)  triangular  plate  element.  A  common  feature  of 
all  of  these  methods,  mixed  formultions,  reduced  integration  and  mode 
decomposition,  is  that  they  can  be  viewed  as  stress  projections.  Stress 
projections  are  methods  in  which  the  stresses  are  projected  on  a  subspace  of 
stresses.  If  the  stress  projection  is  designed  so  that  parasitic  shear 
stresses  are  reduced,  then  shear  locking  is  mitigated.  Similarly,  stress 
projections  that  reduce  parasitic  membrane  stresses  reduce  membrane  locking. 

The  classification  of  projection  methods  as  developed  in  this  contract  is 
summarized  in  Fig.  2.  They  are  classified  as  (1)  implicit  projection  methods, 
such  as  reduced  integration,  and  (2)  explicit  projection  methods,  such  as 
mixed  methods  and  mode  decomposition  methods.  Mode  decomposition  projections 
are  the  most  explicit  of  these  methoas  in  that  algebraic  procedures  are  used 
to  remove  parasitic  stresses.  All  of  the  projection  methods  ameliorate 
locking  by  annihilating  parasitic  shear  and  membrane  stresses. 


The  outline  of  the  report  is  as  follows:  in  Section  1  a  simple  model 
will  be  used  to  show  the  similarity  of  the  causes  of  shear  and  membrane 
locking  and  their  relationship  to  parasitic  shear  and  membrane  stresses. 
Section  2  will  then  describe  how  mode-decomposition  stress  projection  methods 
can  be  used  to  alleviate  shear  and  membrane  locking.  Section  3  will  show  that 
in  mode  decomposition  projection  methods  the  standard  B-matrix  is  projected 
onto  the  interpolation  for  the  stresses.  In  Section  3  it  is  shown  that  the  9- 
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node  Lagrange  element  with  2x2  integration  represents  a  stress  projection 
which  is  free  of  parasitic  shear  and  membrane  stresses,  and  hence  free  of 
membrane  and  shear  locking.  This  element  is  rank  deficient  and  requires 
stabilization;  a  y- stabilization  developed  in  this  contract  is  described  in 
Appendix  A.  In  Appendix  B,  a  formal  equivalence  is  established  between  this 
element  and  an  exactly  integrated  mixed  method. 

In  Section  5  a  challenging  set  of  tests  problems  for  linear  analysis  of 
shells  is  desscribed.  We  have  found  these  problems  to  be  very  decisive  in 
establishing  the  viability  of  elements  and  have  called  it  an  obstacle 
course.  The  performance  of  the  new  elements  developed  here  and  some  older 
elements  on  the  obstacle  course  is  described  in  Section  6.  It  is  concluded 
that  the  9-node  element  with  uniform  2x2  quadrature  and  the  y-stabil ization 
developed  in  this  contract  performs  superbly  on  this  set  of  problems. 
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1.  SHEAR  AND  MEMBRANE  LOCKING  IN  CURVED  BEAMS 


For  purposes  of  examining  the  causes  of  shear  and  membrane  locking,  we 
will  consider  a  curved  beam  described  by  a  one  dimensional  version  of  the 
Marguerre  shallow  shell  equations.  Although  this  is  one  of  the  less  popular 
methods  for  treating  curved  shells  by  finite  elements,  it  should  be  stressed 
that  the  mechanical  behavior  of  elements  described  by  alternative  methods  such 
as  degenerate  shell  theories  and  classical  deep  shell  theories  is  identical  as 
long  as  the  shell  element  is  shallow;  the  convergence  of  shallow  shell 
equations  expressed  in  Cartesian  components  to  deep  shell  results  has  been 
effectively  argued  by  Idelsohn  (1981).  Stolarski  et  al .  (1985)  have  extended 
those  arguments  to  degenerated  shell  theories. 

In  most  practical  applications,  shell  elements  are  quite  shallow  because 
larger  elements  would  prevent  the  achievement  of  satisfactory  accuracy. 
Furthermore,  locking  effects  Increase  with  increasing  curvature.  Thus  the 
Marguerre  theory  provides  an  ideal  vehicle  for  the  study  of  locking  in  curved 
elements. 

The  kinematic  relations  for  the  Marguerre  beam  are  given  by 
e  =  u,  +  w?  w,  (1.1) 

<  -  -*,x  (1.2) 

y  a  w,x  -  *  (1.3) 


Here  u  and  w  are  the  x  and  y  components  of  the  displacement  of  the  midline 
and  $  is  the  rotation  of  the  cross-section  and  the  rigid  body  motion  is 


t 
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removed  from  w,  so  that  w  is  the  displacement  relative  to  the  chord,  see 
Belytschko  and  Glaum  (1979);  w°  is  the  initial  deflection  of  the  midline  from 
the  chord  of  the  element  c,  <  and  y  are  the  membrane  strain,  curvature  and 
transverse  shear  strain,  respectively;  x  is  the  chord  of  the  element  (see  Fig. 
3)  and  commas  denote  differentiation. 

The  stiffness  matrix  for  an  element  is  obtained  in  the  displacement 
method  by 


f  a  3U 

3^ 


(1.4) 


where  dj  is  a  nodal  degree  of  freedom,  f^  the  corresponding  nodal  force  and  U 
the  potential  energy.  For  an  elastic,  isotropic  beam  the  potential  energy  is 
given  by 


U  *  j  f  (Dbk2  +  Dm£2  +  0sy2)  dfl 


(1.5) 


where  Q  is  the  domain  of  the  element  (length  L) ;  Dg ,  DM  and  D«.  are  the 
bending,  membrane  and  shear  constants.  For  an  elastic  beam  of  thickness 
(depth)  d  and  unit  width,  Young's  modulus  E  and  shear  modulus  G,  these 
constants  are  given  by 


»B  *  T7  “3 


(1.6a) 


0M  *  Ed 


(1.6b) 


D 


S 


<sGd 


(1.6c) 


A'  -V  . 
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where  «s  is  the  shear-correction  factor.  The  energies  associated  with  the 
constants  Dg,  and  Dg  are  called  the  bending,  membrane  and  shear  strain 
energies.  Note  that  for  a  thin  beam,  d/L  «  1, 


°B  °B 

«  i  —i.  «  i 

L  °m  L  DS 


(1.7) 


and  the  ratios  of  bending  to  membrane  and  bending  to  shear  energy  will 


similarly  be  small  if  the  strains  are  of  equal  order. 


Shear  and  membrane  locking  arises  in  curved  elements  because  of  the 
inability  of  most  finite  elements  to  achieve  deformed  states  in  which  the 
transverse  shear  strain  and  membrane  strain  vanish  throughout  the  element. 


5 


5 


Modes  of  deformation  in  which  shear  and  membrane  stresses  vanish  play  an 
important  role  in  the  mechanics  of  shells.  For  example,  when  the  moment  field 
is  constant,  the  transverse  shear  vanishes.  When  a  shell  with  a  single 
nonzero  curvature,  i.e.  a  cylindrical  shell,  is  subjected  to  a  state  of  pure 
bending,  the  membrane  strains  will  vanish.  This  mode  of  deformation  is  called 
an  inextensional  mode  of  deformation  because  when  the  membrane  strain 


vanishes,  all  lines  in  the  middle  surface  of  the  shell  remain  constant  in 
length. 

In  a  curved  finite  element,  inextensional  states  of  deformation  are  often 
not  possible.  The  consequences  of  this  shortcoming  are  severe  when  finite 
elements  are  used  to  analyze  a  shell  which  undergoes  inextensional  bending. 
Because  of  the  Inequalities  which  hold  for  thin  beams,  (1.7),  even  a  small 
membrane  strain  or  transverse  shear  strain  will  cause  the  membrane  or  shear 
energy  to  overshadow  the  bending  energy.  Therefore,  if  a  finite  element  is 
used  to  model  a  shell  deforming  in  pure  bending,  it  must  be  capable  of 
representing  this  deformation  so  that  only  the  bending  energy  is  nonzero.  Any 


shear  or  membrane  strains  which  are  developed  will  absorb  a  substantial  amount 
of  energy  and  the  element  will  behave  too  stiffly,  which  is  known  as 
locking.  The  stresses  associated  with  these  spurious  energies  are  often 
called  parasitic  shear  and  membrane  stresses.  Elimination  of  parasitic  shear 
and  membrane  stresses  will  eliminate  locking. 

It  will  now  be  shown  that  many  of  the  commonly  used  elements  will  exhibit 
locking.  Substituting  Eqs.  (1.1)  to  (1.3)  into  (1.5)  yields 

u  3  7  /  [Dr$»v  ♦  bending  energy 

c  a  0 

+  Dm(u,x  +  w’x  W*J2  *■  metnbrane  energy  (1.8) 

+  °s( w*x  ■  4>J  ]  dfl  +  shear  energy 

If  we  consider  a  quadratic.  Isoparametric  element,  then  the  three  displacement 
components  and  the  initial  displacement  of  the  midline  from  the  chord  are 
given  by  quadratic,  Lagrange  shape  functions  ,  so 

3 

[u,  w°,  $,  w]  =  j  Nj(x)  [Uj,  w°,  <t>x ,  Wj]  (1.9) 

From  Eqs.  (1.8)  and  (1.9),  it  can  be  seen  that  if  (1)  w°  *  0  (i.e.  if  the 
shell  is  curved  -  see  Fig.  3),  and  (2)  if  w  *  0,  then  the  membrane  energy  will 
be  nonzero  because  u,x  is  linear  and  cannot  cancel  w°x  w,x,  which  is 
quadratic.  Since  the  two  conditions  of  the  previous  sentence  are  met  in  a 
state  of  inextensional  bending,  this  element  will  exhibit  parasitic  membrane 
stresses  and  be  subject  to  membrane  locking  when  the  stiffness  is  integrated 
accurately;  this  behavior  has  been  studied  by  Stolarski  and  Belytschko  (1983). 
A  9-node  Lagrange  shell  element  will  encounter  the  same  difficulties. 


Remark .  Note  that  w  is  the  transverse  displacement  relative  to  the  chord,  so 
that  it  vanishes  in  rigid  body  motion;  otherwise,  nonzero  membrane  strain 
appears  in  rigid  body  rotation. 

Similarly,  because  $  is  a  polynomial  of  higher  order  than  w,x  in  this 
element,  it  is  clear  from  (1.8)  that  parasitic  shear  energy  will  appear. 
However,  because  the  order  of  the  polynomial  associated  with  parasitic  shear 
is  lower  than  that  associated  with  parasitic  membrane  strains,  shear  locking 
will  not  be  as  pervasive  or  damaging  in  this  element  as  membrane  locking. 

Similar  locking  mechanisms  can  also  be  shown  to  occur  for  the  cubic, 
isoparametric  Lagrange  elements,  although  the  results  of  Arnold  (1981)  show 
that  the  locking  phenomenon  diminishes  as  the  order  of  the  polynomials 
increase.  For  linear  isoparametrics,  no  membrane  locking  occurs  because  w°x 
vanishes. 

In  Kirchhoff  or  C1  elements,  only  membrane  locking  occurs  since  the 
relation  w,x  *  41  is  incorporated  in  the  element,  so  that  the  transverse  shears 
vanish.  The  membrane  locking  phenomenon  in  these  elements  can  be  quite 
severe.  For  example,  in  a  standard  beam  with  w°  and  w  interpolated  by  cubic 
Hermite  interpolants  and  u  interpolated  by  linear  interpolants,  the  term 
w°x  w,x  is  a  quartic,  so  it  cannot  be  effectively  negated  by  u,x,  which  is 
only  a  constant. 

Some  Investigators  have  advocated  using  a  higher  order  interpolation  for 
u  than  w  so  that  the  membrane  strain  can  be  eliminated  in  pure  bending  modes. 
For  example,  for  cubic  w°  and  w,  a  quintic  interpolant  for  u  would  enable 
spurious  membrane  strains  to  be  suppressed.  However,  this  remedy  introduces 
substantial  drawbacks  in  nonlinear  analysis,  for  as  shown  by  Argyris  and 
Scharpf  (1968),  unless  the  shape  functions  are  of  the  same  order  as  that  used 


to  represent  the  geometry,  difficulties  arise  in  large  rotations.  Rigid  body 
motion  can  be  adequately  represented  with  shape  functions  of  different  order 
when  a  corotational  coordinate  system  is  used  for  the  element,  as  in 
Belytschko  and  Hsieh  (1973).  However,  the  use  of  a  single  corotational  system 
for  an  element  with  high  order  shape  functions  may  introduce  errors  because 
rigid  body  rotation  then  varies  substantially  within  a  single  element. 

In  summary,  if  the  orders  of  the  interpolation  polynomials  for  u  and  w 
are  such  that  constant  moment  states  generate  shear  and  membrane  energy,  that 
is  parasitic  stresses,  then  locking  will  result.  Removal  of  the  parasitic 
stresses  is  therefore  a  remedy  for  locking. 

Remark  1.1.  Elements  which  exhibit  locking  do  ultimately  converge.  Locking 
does  not  imply  the  absence  of  convergence,  but  indicates  the  inability  of  the 
element  to  provide  reasonable  accuracy  for  coarse  and  moderate  meshes. 

Remark  1.2.  Many  alternative  paradigms  are  available  for  locking.  For 
example,  Prathap  and  Bhashyam  (1982)  explain  locking  by  the  appearance  of 
spurious  constraint  equations;  Park  and  Flaggs  (1984)  explain  locking  through 
the  appearance  of  unusually  high  frequencies  in  a  Fourier  analysis  of  the 
discrete  equations.  All  of  these  approaches  provide  useful  insights  into  the 
causes  of  locking,  but  in  the  setting  of  projection  methods,  the  viewpoint  of 
parasitic  stresses  appears  most  elucidating. 


2.  MODE  DECOMPOSITION  PROJECTION  METHODS 


The  purpose  of  the  mode-decomposition  projection  methods  is  to  project 
the  nodal  displacements  so  that  the  parasitic  stresses,  and  hence  locking,  are 
eliminated.  In  this  Section,  the  application  of  projections  to  eliminate 
shear  and  membrane  locking  in  simple  beam  elements  will  be  described;  the 
application  of  projection  to  a  triangular  shell  element  is  also  sketched. 

The  basic  idea  of  mode  decomposition  projection  methods  is  to  define  the 
bending  mode  component  of  any  deformation  and  to  ignore  the  membrane  and  shear 
strain  energies  associated  with  the  bending  modes.  Thus  the  strain  energy 
expression  (1.5)  becomes 

U  *  j  f  [Db»c2  ♦  0M(e  -  cb)2  *  Ds(y  -  Yb)2]dQ  (2.1) 

where  eb  and  yb  are  the  membrane  and  shear  strains  in  the  bending  mode. 

As  can  readily  be  seen  from  the  above,  in  a  constant  moment  state,  the  bending 
mode  constitutes  the  total  deformation  so  the  membrane  and  shear  energies  will 
vanish,  since  in  that  case  e  *  eb  and  y  *  yb. 

This  mode  decomposition  is  implemented  by  defining  bending  nodal 
displacements  through  a  projection 


The  requirement  that  the  bending  parts  of  the  membrane  and  shear  strains 
£b  and  rb,  be  those  strains  that  occur  in  a  constant  moment  state  is  not 
sufficient  to  identify  the  projection  operator  Pfa,  because  it  is  only  one  of 
the  bending  modes.  For  other  nodal  displacements,  it  is  also  necessary  to 
define  the  bending  mode  (or  component)  of  the  nodal  displacements.  Therefore 
a  more  general  procedure  must  be  developed. 

Let  us  consider  the  mode  decomposition  procedure  for  shear  projection  in 
C°  beam  elements.  The  first  step  of  this  procedure  is  to  define  a  Kirchhoff 
mode  for  any  displacement  field.  A  Kirchhoff  mode  is  the  displacement 
field  w  (x)  that  would  occur  in  a  Kirchhoff  beam  with  the  curvature  field 
given  by  Eq.  (1.2);  w  (x)  is  obtained  by  simply  Integrating  the  curvature 
twice  and  then  choosing  the  constant  of  integration  so  that  the  Kirchhoff 
mode  best  fits  the  nodal  displacements  of  the  element. 

For  example,  consider  the  2  node,  linear  w,  linear  $  element  in  which 


W  *  WjU  -  5)  +  w2  £ 

(2.4a) 

*  ■  ^(1  -  5)  +  <>2  5 

(2.4b) 

In  this  element,  Eqs.  (1.2)  and  (2.4b)  yield  the  following  curvature 

*  *  ■  r  "  *1^  1  ”  *21^  (2.5) 

and  the  Kirchhoff  modes  are  given  by 


The  constant  is  selected  by 


-  .  (2.7) 

so  that  the  nodal  rotations  of  the  Kirchhoff  mode  match  the  total  nodal 
rotations  of  the  beam 


w’x(0)  *.*1 

(2.8a) 

w,xK(L)  -  *2 

(2.8b) 

This  selection  is  based  entirely  on  physical  reasoning;  the  aim  is  to  insure 
that  in  a  constant  moment  state,  such  as  shovwi  in  Fig.  3,  the  rotations  are 
entirely  associated  with  the  Kirchhoff  (i.e.  the  bending)  mode,  so  that 
transverse  shear  energy  vanishes  in  (2.1).  It  is  not  always  possible  to 
satisfy  relations  such  as  (2.8)  (see  for  example  the  development  of  the 
triangular  plate  in  Belytschko,  et  al .  (1984)). 

The  Kirchhoff  mode  is  then  selected  as  the  bending  mode,  so  that  the 
relation  between  the  total  nodal  displacements  and  the  bending  nodal 
displacements  follows  from  (2.6-8)  and  is  given  by 


A  little  manipulation  of  Eqs.  (2.1),  (2.3)  and  (1.4)  shows  that  the  stiffness 
of  this  element  is  given  by 


where 


where  £  is  a  unit  matrix;  Eqs.  (1.1-3)  and  (2.4)  give 

Bb  *  |  [-1,  +1,  0,  0] 

Bs  ■  [5  -  1,  -e,  -  l/l,  1/L] 


(2.10a) 

(2.10b) 

(2.11a) 

(2.11b) 


An  Interesting  consequence  of  the  fact  that  £s  is  a  projection  operator  is 
that  (see  Stolarski  and  Belytschko  (1985)) 


(2.13a) 


K*  -  Pi  /  b:  Oe  8  dn  Pc 
~e  ~s  ‘  ~s  ~s  ~s  ~s 


?!  /  BL  De  B  dn  Pc 
~S  jj  ~ss  ~s  ~ss 

e 


C  bt  0  8  dn 

1  ~ss  ~s  ~ss 


where 


(2.13b) 

«.*\*\* 

>V, 

■>v 

V-.v 

f. 

(2.13c) 

Sss  3  Ss£s  3  t-  1/2’  -1/2*  -1/L‘  +1/L^ 


(2.14) 


The  effect  of  the  projection  is  thus  to  change  the  shear  strain  from  a  linear 
field  such  as  given  by  the  Bs-matrix  of  Eq.  (2.11b)  to  the  constant  field 
given  by  (2.14) . 

Remark  2.1.  It  can  be  seen  from  (2.13c)  and  (2.14)  that  the  stiffness  matrix 
obtained  by  this  procedure  is  identical  to  that  obtained  by  reduced/selective 
integration  by  Hughes,  et  al .  (1977). 

Remark  2.2.  Considerable  leeway  is  available  in  the  choice  of  the  bending 
mode.  We  have  chosen  to  ascribe  all  of  the  nodal  rotations  to  the  bending 
mode,  but  other  choices  are  acceptable. 

Remark  2.3.  The  linear  displacement,  linear  rotation  element  is  not  subject 
to  membrane  locking,  since  w°x  w,x  vanishes  in  (1.8)  for  that  element. 

The  projection  method  for  membrane  locking  is  similar.  In  this  case,  the 
aim  is  to  find  the  membrane  strain  eb  so  that  in  an  inextensional  mode  of 
deformation,  the  membrane  energy  vanishes  in  Eq.  (2.1).  The  projection  for 
the  linear-cubic  beam  element  has  been  developed  by  Stolarski,  et  al .  (1984a), 


(1984b)  and  related  to  a  -nixed  method  by  Selytschko,  et  al  .  (1985b).  The 
displacement  a  o  initial  deflection  in  this  element  are  given  by 


u  »  u.(l  -  C)  +  u0  Z 


(2.15a) 


w  -  ♦j  L  (C3  -  252  +  C)  +  *2  L  (53  -  52) 


(2.15b) 


»?  L(C3  -  2?2  +  O  +  L(C3  -  C2) 

*  V  J  £  v  7 


(2.15c) 


Note  again  that  we  have  used  a  corotational  coordinate  system  so  that  the  x 
axis  always  connects  node  1  to  node  2. 

The  stiffness  matrix  for  this  element  is  obtained  by  using  (1.1-2)  in 
conjunction  with  (2.1)  which  yields 


K  *  /  3J  Dn  B.  dn  +  ?l  /  B  Dm  Bm  dn  P 

~e  ^  ~b  ~B  ~b  ~m  '  — m  ~m  ~n 

e 


(2.16) 


where 


Sb  *  tO. 


(2.17a) 


tin  “  It  ’  "°(x)Nl,x-  w“(x)"2,xl 


(2.17b) 


d  a  [u£  ■  u^,  4>  ^ »  ^2! 


(2.17c) 


If  =  l,  or  in  other  words,  if  no  projection  is  used,  then  this  element  will 
lock.  The  cause  of  locking  can  be  understood  by  simply  considering  3  ;  if  the 


moment  on  an  element  is  constant,  then  t> ^  so  if  w  (x)  *  0,  then  the 

constant  term  1/L  cannot  negate  the  terms  w°(x)N.  ,  hence  parasitic  membrane 

1  ,  X 

strains  will  arise.  If  w°(x)  =  0,  the  beam  is  straight  and  no  coupling  exists 
between  flexural  and  membrane  effects.  This  coupling  is  an  important 
attribute  of  curved  elements. 

The  projection  operator  which  eliminates  locking  is  developed  in 
Stolarski,  et  al .  (1984)  to  be 

1  Iff  (_4*l  +  ffff  <♦?  "  ^ 

00  0  '  (2.18) 

0  0  0 

This  projection  is  obtained  by  noting  that  the  change  in  the  chord  length 
U21  S  u2  "  U1  ln  a  pure  bendln9  11,0(16  ls  given  by 

U21  =  J  '  W°  W,xx  dx  =  Iff  (-4*l  +  *2^1  +  ffff  (*1  ‘  4*2^2  (2*19) 

As  can  be  seen  from  Fig.  3,  this  change  in  chord  length  is  required  to 
maintain  an  inextensible  midline  during  bending  of  a  curved  element;  if  an 
element  is  straight,  no  change  in  chord  length  occurs  during  bending. 

An  analogous  approach  has  also  been  used  to  formulate  a  triangular  shell 
element  In  Stolarski,  et  al .  (1984)  which  is  here  called  DKT-CST*.  In  its 
development,  it  is  convenient  to  use  only  deformational  degrees  of  freedom. 
Following  Argyris  (1965),  the  membrane  state  of  strain  is  expressed  in  terms 
of  the  elongations  of  the  3  sides  n^  and  the  bending  state  represented  by 
deformation  rotations  VlT,  so  that  the  deformational  d.o.f.  are 


i 

*•? 


.W"  Jtt  J 


fc.ffc.4- 1 


il.ililU 


nT  -  [nx,  n2,  n3] 


i  3  C*xl  •  *yl’  *x2>  V’  ^x3  * 


The  strain  energy  of  the  element  is  then  given  by 


u  3  \  &T  JSb  &  +  7  Cl!  '  nh)T  <m(n  -  nb) 


(2.20a) 


(2.20b) 


b'  ~m' 


(2.21) 


where  Kb  and  <m  are  the  bending  and  membrane  stiffnesses,  respectively.  On 

the  element  labeled  DKT-CST*,  in  Section  6,  the  discrete  Kirchhoff  triangle 

* 

(DKT)  element  of  Batoz  (1982)  is  used  for  the  bending  stiffness  and  the 
constant  strain  element  is  used  for  membrane  stiffness.  Note  that  in  the 
absence  of  ^  in  (2.21),  the  bending  and  membrane  behavior  of  the  element  are 
totally  uncoupled. 

The  bending  elongation  along  each  side  is  now  given  by  the  same  equations 
as  for  the  beam,  (2.19),  but  written  in  the  form 


nIb  =  17  II  +  *2l)  +  ll  +  17  ^11  '  4^2I^2 


(2.22) 


where  and  are  the  rotations  relative  to  side  I  of  the  element  and 

*11’  ^21  are  ttle  slopes  of  the  shell  surface  relative  to  the  chord  of 

the  element.  From  the  relations  (2.22)  at  the  three  sides,  the  projection 

operator  P_  is  obtained. 

~m 

Details  of  the  element  are  given  in  Stolarski,  et  al ;  (1984).  One 
important  aspect  of  the  formulation  which  was  inadvertently  omitted  from  that 
paper  is  that  and  <p^  in  (2.22)  must  be  defined  so  that  the  bending 


elongations  of  the  common  edge  of  two  contiguous  elements  are  identical. 
Otherwise,  the  projections  for  elements  which  are  not  coplanar  will  be 
different,  and  as  a  consequence,  parasitic  membrane  energy  will  appear  in  a 
state  of  pure  bending.  This  is  accomplished  by  defining  ^  as  the  rotation 
in  the  plane  defined  by  side  I  and  the  normal  to  the  shell  at  node  i. 

It  can  be  seen  from  (2.22)  that  the  bending  part  of  the  elongation 
couples  the  membrane  response  with  the  rotations,  and  hence,  as  in  the  beam, 
adds  membrane/flexural  coupling  to  the  element. 


3.  RELATION  BETWEEN  MIXED  AND  MODE-DECOMPOSITION  PROJECTION  METHODS 


In  this  Section  a  simple  form  of  the  element  stiffness  matrix  for  mixed 
methods  will  be  developed  and  then  compared  to  the  explicit  projection  method 
to  illustrate  their  similar  structure.  The  Hu-Washizu  variational  principle 
will  be  used  to  develop  the  mixed  method.  This  principle  can  be  written  for  a 
single  element  in  the  form 


[  eu  -  «ij)  -  s°ijhj  - 

nN 

*  5uo,j)  oijld0  *  £  iun  fu 


(3.1) 


Standard  indiclal  notation  has  been  used,  with  repeated  lower  case  subscripts 
summed  over  their  range;  commas  denote  partial  derivatives.  Upper  case 
subscripts  pertain  to  degrees  of  freedom  of  the  nN  nodes.  The  nomenclature  in 
this  section  is  ,as  follows: 

e..  ,  e  =  strain;  matrix  form  e  is  considered  a  column 

1 J  *** 

matrix  such  as  (ex,  ey,  2exy) 

a . .  ,  a  «  stress  matri x 

ij 

u.  ,  u  *  displacement  field 

dj  ,  d  *'  nodal  displacements 

u^  ,  vsu  =»  symmetric  part  of  the  gradient  of  the 

displacement  field 

^ i j ki  *  2  *  ■  stress-strain  matrix 

In  the  above,  both  the  subscripted  tensor  forms  and  the  matrix  forms  commonly 
used  in  f i n i element  implementation  are  listed.  Through  the  use  of  the 


latter  form,  the  relation  of  the  projection  method  and  Hu-Washizu  can  be 
clarified  more  easily.  The  displacement  field  in  this  principle  must  be  C°, 
while  the  strain  and  stress  fields  may  be  C"*. 

The  three  independent  fields  are  approximated  by  interpolation  functions 
as  follows 

nD 

u4  3  T  N . T  dT  or  u  3  N  d  (3.2a) 
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where  nQ  are  the  number  of  displacement  d.o.f.  and  n^  are  the  number  of  strai 
(or  stress)  interpolants.  An  important  part  of  this  presentation  is  the  use 
of  mnemonic  terms  for  the  strain  and  stress  interpolants  and  coefficients  so 
that  they  are  easily  recognized  later.  We  also  define  the 
standard  B-matrix  through 
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Substituting  the  above  into  Eq.  (3.1),  and  using  the  arbitrariness  of  the 
appropriate  variations  yields: 


strain-displacement  equations 


nE  _  n0  _ 

T  E. .  e  ,  3  l  B. ,  d  , 
IJ  J  ,  ,  IJ  J 
J*1  J*1 


or 


E  e  3 


stress-strain  equations 
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At  this  point  it  is  convenient  to  note  that 
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In  terms  of  this  notation,  (3.6)  can  be  rewritten  as 
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From  this  notation  it  becomes  clear  that  the  essential  ingredients  of  the 
mixed  method,  such  as  the  8-matrix,  are  simply  the  projections  of  these 
matrices  onto  the  corresponding  stress  and  strain  interpolants.  The 
^-matrix,  for  example,  is  projected  onto  the  stress  interpolants  in  this  mixed 
method. 

A  more  revealing  form  of  the  mixed  method  can  be  obtained  by 
orthogonal  1  zing  the  strain  and  stress  interpolants  so  that 
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where  _I  is  the  unit  matrix.  This  can  always  be  achieved  by  a  Gram-Schmidt 
procedure  if  the  stress  and  strain  interpolants  span  the  same  space  and  are 
linearly  independent;  it  is  also  convenient  if  each  strain  (or  stress) 
parameter  (or  s^)  pertains  only  to  a  single  strain  (or  stress  o^). 
Using  (3.9),  Eqs.  (3.3  -  5)  can  be  written 
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which  can  easily  be  recognized  as  the  strain  -  displacement,  constitutive  and 
nodal-force  stress  (equilibrium)  relations  in  projected  form. 

The  element  stiffness  matrix  is  now  obtained  by  combining  (3.4-6),  which 
is  most  easily  done  with  the  matrix  form,  and  gives 

f  *  BT  (E1)"1  D  (E"1)  B  d  (3.11) 
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Using  the  orthogonality  (3.9),  (3.11)  can  then  be  written  as 

Se'Ffil  (3.12) 

For  purposes  of  comparison  of  the  mixed  method  with  the  mode  decomposition 
projection  method  developed  in  Section  2,  it  is  also  worthwhile  to  write 
(3.12)  using  the  definition  of  Ojj,  (3.8c),  which  gives 

Ko  -  BT  /  ET  0  E  dn  *8  (3.13) 

a 

Table  1  compares  the  above  form  with  the  mode  decomposition  projection 
method.  It  can  be  seen  that  the  two  methods  are  identical  in  structure,  with 
the  matrix  P  playing  the  role  of  the  projection  of  g  onto  the  stress 
interpolant  S,  hence  the  name  stress-projection . 

In  many  elements,  the  same  interpolation  functions  are  used  for  the 
stresses  and  strains.  To  establish  the  equivalence  of  the  two  methods,  the 
strain  interpolants  should  be  orthogonal i zed ,  and  the  matrix  P  plays  the  role 
of  the  projection  of  the  B  matrix  onto  the  strain  interpolant  r. 


L  f . 


Although  the  conditions  for  equivalence  between  mode  decomposition 
projection  methods  and  mixed  methods  appears  straightforward,  establishing  the 
exact  equivalence  is  not  possible  in  some  elements.  Stolarski  and  Belytschko 
(1985a)  have  shovwi  an  equivalence  between  mode  decomposition  projection 
methods  and  mixed  methods  for  a  large  variety  of  beam  elements,  but  in  that 
paper  it  was  shown  that  no  mixed  element  is  equivalent  to  the  triangular  plate 
element  developed  in  8elytschko,  et  al .  (1984c). 

Remark .  The  equivalence  of  the  two  methods  does  not  depend  on  the 
orthogonal ity  of  the  strain  and  stress  interpolants ;  it  only  serves  to  clarify 
the  relationship.  The  stiffness  matrix  in  fact  is  not  changed  by 
orthogonal i zati on  of  the  interpolants. 

Table  1 

Comparison  of  Mode  Decomposition  Projection  and  Mixed  Methods 

Mixed  Method 

(Orthogonal  Interpolants) 

Strains 

See  Eqs .  (3.2b),  (3.10a) 

£  *  £  S.  -  l  a  <t 

Stiffness  Matrix 

See  Eq.  (3.13)  See  Eq.  (2.10a)  (only  shear  stiffness) 

k  s  3T  /  ET  0  E  dn  3  k  -  pT  /  b!  D„  da 

Note  interchangeabi 1 ity  when  8  =  P  ,  E  =  Bs . 


4.  NINE-NODE  LAGRANGE  SHELL  ELEMENT 


The  9-node  shell  element  to  be  discussed  here  has  been  extensively 
studied  by  Parish  (1979)  and  by  Hughes  and  Liu  (1981).  It  is  a  Mind! in  type, 
degenerated  shell  element  which  uses  quadratic  Lagrange  interpolants.  The  2x2 
quadrature  version  of  the  9-node  element  possesses  the  unique  and  beneficial 
property  that  in  constant  moment  states,  the  transverse  shear  and  membrane 
stresses  (that  is,  the  parasitic  stresses),  vanish  at  the  quadrature  points. 
Thus  2x2  quadrature  provides  a  stress  projection  which  should  avoid  locking. 

This  property  can  be  verified  by  performing  the  numerical  experiments 
illustrated  in  Fig.  4.  As  indicated  in  the  figure,  both  flat  elements  and 
elements  with  a  single  curvature  were ‘considered .  In  the  numerical 
experiment,  moments  were  applied  as  shown  and  the  transverse  shear  energy  and 
membrane  energy  was  monitored.  In  all  cases  which  were  tried,  the  transverse 
shear  and  membrane  energies  were  less  than  0.01%  of  the  total  strain  energy. 

The  same  behavior  was  observed  when  the  quadrature  points  were  shifted 
from  n  =  5  3  t  3“  ^  (the  Gauss  points)  to  n  *  £  =  i  V2  •  However,  the  accuracy 
of  the  element  then  deteriorates  even  though  the  membrane  energy  locking  does 
not  occur;  for  an  illustration  of  this  see  Table  2,  which  gives  results  for 
the  arch  problem  described  in  Fig.  lb.  These  runs  were  made  with  a  beam 
element  because  the  shell  develops  a  spurious  mode  for  the  ±  V2  quadrature 
scheme. 

This  feature  of  the  9-node  element  with  2x2  quadrature  makes  it  a  very 
desirable  element  for  the  analysis  of  shells.  Its  highly  convergent  behavior 
will  be  illustrated  in  Section  6.  However,  the  element  suffers  from  one 
important  drawback:  its  rank  is  not  sufficient  to  preclude  spurious  singular 
modes,  so  for  some  boundary  conditions  the  total  stiffness  matrix  is  singular. 
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Table  2 

Effect  of  Quadrature  Point  Location  on  Accuracy  for  the  Arch  Problem 

Number  of  Elements  wFEM/wexact  U  /U  U  /U 

m  s 


.l/o 

Gauss  quadrature  5  *  ±  3  ‘ 


10 

0.9861 

0.0013 

0.0031 

20 

0.9984 

0.0013 

0.0031 

40 

1.0009 

0.0013 

0.0031 

Midpoint  quadrature  5  =  ±  1/2 


10 

1.0928 

0.0016 

0.8033 

20 

1.0257 

0.0120 

0.0041 

40 

1.0077 

0.4199 

0.0031 

Um  3  membrane  energy,  Us  3  shear  energy,  U  3  total  energy 


VvV 


s*\.* 

V<* 


w  3  deflection  under  load 


•The  control  of  these  spurious  modes  has  been  developed  for  the  9-node 
plate  element  in  Selytschko,  et  al .  (1984),  and  the  stabi 1 i zation  procedure 
for  the  shell  element  has  been  reported  without  analysis  in  Belytschko,  et  al 
(1985).  In  Appendix  A  a  description  of  the  element  as  used  here  and  a 
detailed  development  of  the  stabilization  operator  is  given. 

Reduced  integration  of  the  stiffness  may  also  be  viewed  as  a  projection 
method.  This,  of  course,  is  suggested  by  the  equivalence  principle  of  Malkus 
and  Hughes  (1978)  which  establishes  the  equivalence  of  reduced  integration 
displacement  and  reduced  integration  mixed  methods;  the  relationship  between 
the  latter  and  projection  methods  has  already  been  discussed  herein.  Reduced 
integration  represents  an  implicit  projection,  in  that  the  projection 
operator  is  never  explicitly  invoked.  An  exact  mixed  formulation  for  the  9- 
node  (2x2)  element  is  developed  in  Appendix  B.  In  that  development,  by  using 
Dirac  functions  for  the  stresses,  the  reduced  quadrature  displacement 
formulation  can  be  shown  to  be  equivalent  to  an  exactly  integrated  mixed 
element. 


5.  AN  08STACLE  COURSE  FOR  SHELL  ELEMENTS 


A  major  shortcoming  which  has  impaired  the  development  of  shell  elements 
is  that  in  the  reports  on  many  elements,  they  are  not  tested  against  a  set  of 
problems  which  challenge  even  a  fraction  of  the  capabilities  required  in  a 
high-performance  element.  A  good  shell  element  must  have  the  ability  to 
handle  inextensional  bending  modes  of  deformation,  rigid  body  motion  without 
straining,  and  complex  membrane  states  of  stress.  Inadequacies  in  this 
spectrum  of  attributes  are  a  severe  handicap. 

A  useful  obstacle  course  for  an  element  must  also  be  reasonably  short;  it 
is  useless  to  include  problems  which  only  disqualify  elements  which  are 
already  disqualified  by  other  problems. 

In  this  effort,  we  have  assembled  the  three  test  problems  shown  in  Fig. 

5.  All  of  these  problems  have  been  seperately  used  by  others  in  the 
literature  as  noted  in  Table  3.  We  have  found  that  together  they  are  an 
extremely  discriminating  set  of  problems. 

Some  remarks  on  this  selection: 

1.  The  Scordelis-Lo  problem  is  extremely  useful  for  determining  the  ability 
of  an  element  to  accurately  solve  complex  states  of  membrane  strain.  A 
substantial  part  of  the  strain  energy  is  membrane  strain  energy,  so  the 
representation  of  inextensional  modes  is  not  crucial  in  this  problem. 

Even  elements  with  severe  membrane  locking  will  converge  at  a  moderate 
rate  in  this  test,  whereas  inadequacie-'  in  membrane  stress  accuracy  will 
severely  inhibit  convergence. 

2.  The  pinched  cylinder  with  a  diaphragm  is  one  of  the  most  severe  tests  for 
both  inextensional  bending  modes  and  complex  membrane  states.  We  have  not 


included  the  pinched  cylinder  with  free  ends  because  any  element  that 
passes  the  diaphragm  support  test  will  perform  well  when  the  boundary 
condition  is  simplified  to  a  free  boundary. 

3.  The  hemispherical  shell  problem  is  again  a  challenging  test  of  an 

element's  ability  to  represent  inextensional  modes;  it  exhibits  almost  no 
membrane  strains.  The  role  of  this  test  problem  is  less  critical  with 
regard  to  inextensional  bending  than  the  pinched  cylinder  problem. 
However,  it  is  a  very  useful  problem  for  checking  the  ability  of  the 
element  to  handle  rigid  body  rotations  about  normals  to  the  shells 
surface.  Large  sections  of  this  shell  rotate  almost  as  rigid  bodies  in 
response  to  this  load,  so  that  the  ability  to  accurately  model  rigid  body 
motion  is  essential  for  good  performance  in  this  problem.  Some  5  d.o.f. 
per  node  formulations  of  triangular  elements  fail  this  test  because  they 
result  in  spurious  straining  when  rotated  about  the  normal  to  the  shell 
surface.  This  problem  is  much  more  challenging  than  the  point-loaded 
spherical  problem  which  is  often  used. 


Table  3 


Problem  Parameters  for  Obstacle  Course 


Problem  1.  Scordelis  -  Lo  Roof 

length:  L  *  50.0 
radius:  R  *  25.0 
thickness:  t  *  0.25  a 

Young's  modulus:  E  *  4.32  *  10s 
Poisson’s  ratio:  v  =  0.0 

boundary  conditions:  supported  at  each  end  by  rigid 
diaphragms 

loading:  uniform  vertical  gravity  load  of  90.0  per  unit 
area 

converged  numerical  solution:  vertical  displacement  at 
midside  of  free  edge  =  0.3024 
reference:  Scordelis  and  Lo  (1969),  Ashwell  (1976) 


Problem  2. 


Problem  3. 


Pinched  Cylinder  with  Diaphragms 

length:  L  =  600.0 

radius:  R  =  300.0 

thickness:  t  =  3.0 

Young's  modulus:  E  *  3.0  *  10° 

Poisson's  ratio:  v  =  0.3 

boundary  conditions:  constrained  at  each  end  by  a 
rigid  diaphragm, u  =  u  =  <j>  =  0  in  Fig.  5 
loading:  opposing  radial  ioadszas  shown  in  Fig.  5, 

F  •  UO 

radial  displacement  at  point  load:  0.18248  *  10“4 
references^  lindberg,  et  al .  (1969),  Dvorkin  and  Bathe 
(1984),  FI  ugge  (1973) 

Hemispherical  Shell 

radius:  R  -  10.0 

thickness:  t  *  0.04 

Young's  modulus:  E  3  6.825  *  107 

Poisson's  ratio:  v  =  0.3 

boundary  condition:  bottom  circumferential  edge  of 
hemisphere  is  free 

loading:  opposing  radial  poi’nt  loads  alternating  at  90° 
as  shown  in  Fig.  5,  F  *  ±2.0 
solution:  radial  displacement  at  loaded  points:  0.0924 
references:  Morley  and  Morris  (1978),  MacNeal  and  Harder 
(1984) 


Table  4 


Comparison  of  Results  for  Hemispherical  Shell  (See  Fig.  2) 
with  Consistent  and  Inconsistent  Spurious  Mode  Control 
(results  given  are  ratio  of  computed  to  analytic  displacement  under  load) 


Mesh  for 
Quarter  of  Shell 
(number  of 
nodes/edge) 


Consistent 
Control  y 

given  by  Eq.(A.32) 


1.319 


Inconsistent 
Control  j  *  h 


1.2672 


2x2  Quadrature 

1  s  9 


1.0794 

1.0056 


0.9963 


0.9925 


1.0792 

1.0063 

0.9970 

0.9958 


1.0954 
1.0138 
1  .0037 
0.9987 


Remark  5.1.  This  selection  of  shell  problems  was  inspired  by  the  work  of 
MacNeal  and  Harder  (1984),  who  devised  a  set  of  standard  finite  element 
problems  which  included  hemispherical  shell  and  the  Scordelis  Lo  roof,  but  not 
the  pinched  cylinder.  They  also  included  various  plane  patch  tests  and  plane, 
two  dimensional  problems. 


6.  RESULTS 


Results  have  been  obtained  with  two  classes  of  elements:  3-node 
triangular  elements  with  5  or  6  d.o.f.  per  node  and  C°  Mind! in  type  shell 
elements  based  on  a  degenerate  shell  theory  of  Hughes  and  Liu  (1981).  All  of 
the  triangular  elements  utilize  the  discrete  Kirchhoff  triangle  (DKT) 
formulation  of  Batoz  et  al .  (1982)  with  different  membrane  strain  fields. 

They  are  as  follows: 

1.  DKT  -  CST:  a  standard  flat  triangle  with  constant  membrane  strains  and 
no  membrane/bending  coupling. 

2.  DKT  -  CST*:  a  constant  strain  triangle  is  used  for  the  membrane  strain  as 
above  but  the  membrane  projection  described  in  Section  2  and  originally 
reported  by  Stolarski  et  al .  (1984),  which  couples  membrane  and  bending, 
is  added. 

3.  DKT  -  LST:  a  6  d.  o.  f.  per  node  element  in  which  a  linear  membrane  field 
with  reduced  Integration  as  given  by  Carpenter  et  al .  (1985)  is  employed 
for  the  membrane  effects.  This  element  is  flat  and  has  no  membrane- 
flexural  coupling. 

4.  DKT  -  01 son-Bearden :  a  linear  field  given  by  Olson  and  Bearden  (1979)  is 
used  for  the  membrane  strains  in  combination  with  a  DKT  element.  This 
element  has  6  d.o.f  per  node  and  no  membrane  flexural  coupling. 

Among  the  Mindlin  C°  elements,  the  following  were  used: 

1.  4-node  SRI:  this  is  a  standard  4-node  Mindlin  element  described  by  Hughes 
and  Liu  (1981)  with  selective  reduced  integration,  consisting  of  reduced 
integration  on  shear  (1  point)  and  2x2  quadrature  on  the  bending  and 


membrane  terms. 


2.  9-node  SRI:  the  9-node  Lagrange  element  with  3x3  quadrature  on  all 
bending  and  membrane  terms  and  2x2  quadrature  on  the  shear  terms. 

3.  9-node  3x3:  full  quadrature  (3x3)  on  all  stiffness  terms  of  the  9-node 
element. 

4.  9-node  Y-method:  2x2  quadrature  on  all  terms  of  the  element  stiffness 
with  y-stabil ization ,  as  described  in  Appendix  A. 

Scordelis-Lo  roof.  Results  for  this  problem  shown  in  Fig.  6.  These  are 
observations  of  element  performance  for  this  problem: 

1.  The  DKT-CST  and  DKT-CST*  are  both  very  poor  for  this  problem,  whereas  for 
the  other  problems  in  this  obstacle  course  they  perform  very  well. 

2.  Among  the  triangles,  only  the  01 son-Bearden  and  DKT-LST  elements  perform 
reason  ably. 

3.  The  9-node  Y-element  and  the  4-node  SRI  elements  perform  extremely  well  on 
this  problem.  However,  when  a  selective  reduced  integration  is  used  in 
this  element,  the  results  hardly  differ  from  full  integration.  This  is 
very  puzzling,  since  the  problem  is  dominated  by  membrane  response  rather 
than  bending,  so  the  severe  locking  of  9-node  SRI  element  is  quite 
unexpected.  In  fact,  there  seems  to  be  considerable  membrane  locking  in 
the  coarse  meshes,  but  it  is  quickly  eliminated  with  mesh  refinement. 

Pinched  cylinder  with  diaphragm.  Results  for  this  problem  are  given  in 
Fig.  7.  In  this  problem,  all  of  the  triangular  elements  work  reasonably  well, 
with  only  the  DKT-CST  being  marginal.  On  the  other  hand,  the  performance  of 
all  of  the  Mind! in  elements,  except  the  9-node  Y-element,  are  quite  poor. 
Particularly  noteworthy  is  the  fact  that  even  with  17  nodes  along  each  edge, 
the  4-node  SRI  element  is  still  5%  in  error  and  converging  very  slowly.  This 
is  the  only  problem  in  this  set  in  which  the  4-node  element  performs  below 


Pinched  Cylinder  with  Diaphragm 


par.  No  improvement  is  achieved  in  the  4-node  element  by  using  1-point 
quadrature  on  all  terms,  see  Liu  et  al .  (1985). 

In  order  to  illustrate  the  relative  magnitudes  of  shear  and  membrane 
locking.  Fig.  8  shows  the  results  obtained  for  the  9-node  element  for  SRI, 

2x2,  and  3x3  quadrature.  It  can  be  seen  that  for  coarse  meshes,  membrane 
locking  is  dominant  whereas  for  the  finest  meshes,  shear  locking  becomes  more 
important. 

Figure  9  clarifies  the  role  of  membrane  and  shear  locking  by  showing  the 
fractions  of  the  total  energy  which  are  shear  and  membrane  energies  for 
various  meshes  of  9  and  4-node  elements.  As  can  be  seen,  for  all  of  the 
elements,  the  shear  energy  tends  to  almost  zero  as  the  meshes  are  refined,  but 
the  rate  of  convergence  of  the  shear  energy  for  the  9-node  3x3  element  is 
quite  slow.  The  membrane  energy  fraction  of  the  y-element  tends  to 
approximately  0.373.  For  the  4-node  SRI  element  the  membrane  energy  is 
somewhat  larger  for  all  of  the  meshes  considered. 

For  the  9-node  3x3. element,  the  membrane  energy  is  overpredicted  even  for 
the  most  refined  mesh,  which  is  indicative  of  membrane  locking.  The  initial 
increase  of  the  membrane-energy  fraction,  which  is  apparent  in  Fig.  9,  is 
somewhat  puzzling.  Evidently  for  the  coarse  meshes,  the  total  internal  energy 
is  very  small  because  of  severe  locking;  the  initial  refinement  of  the  mesh  in 
this  element  serves  primarily  to  reduce  shear  locking,  so  that  parasitic 
membrane  energy  actually  increases  at  first. 

Spherical  shell  problem.  The  results  for  the  spherical  shell  are  shown  in 
Fig.  10.  The  following  points  are  noteworthy: 

1.  This  is  the  only  problem  in  which  the  01 son-8earden  element  performs 
poorly. 
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9  node  SRI 
9  node  3x3 
9  node  2x2 


strain  energy  into  membrane  (M)  and  shear  energy  (S)  for  the  9-node  element 
ed  cylinder/diaphragm  and  hemispherical  shell  problems. 


Hemispherical  Shell 


V 


WWWW '.WWW  A 


■ia 


47 


2.  The  other  triangular  elements  'perform  extremely  well  in  this  problem, 
whereas  the  fully  integrated  C°  elements  (4-node  and  9-node)  all  exhibit 
severe  membrane  locking.  Shear  locking  is  almost  totally  absent. 

Remark  6.1:  Noor  and  Peters  (1981)  have  shov*i  similarly  stiff  behavior  of 
displacement  elements  for  curved  beams.  Interesting  results  have  also  been 
given  by  Ramm  and  Stegmuller  (1982)  who  found  that  for  the  16  node  Lagrange 
with  4x4  quadrature,  the  buckling  value  of  a  cylindrical  panel  converged  only 
to  within  2%  of  the  solution  with  a  12x12  element  (37x37  nodes,  6900  d.o.f.) 
mesh;  they  did  not  report  their  curved  9-node  element  results  because  they 
were  so  poor. 

Summary  of  performance.  As  can  be  seen,  this  obstacle  course  includes 
problems  which  compromise  the  performance  of  every  element  except  the  9- 
node  y-element  and  the  DKT-LST.  The  OKT-LST  unfortunately  suffers  the 
drawback  that  it  does  not  have  a  variational  basis.  The  4-node  SRI  element  is 
the  next  best  among  these  elements,  although  it  exhibits  some  locking  in  the 
pinched  cylinder/diaphragm  problem. 

A  confusing  result.  The  senior  author  has  devoted  considerable  care  to  the 
development  of  methods  for  controlling  spurious  modes  which  are  consistent  in 
the  sense  that  the  strains  for  linear  displacement  fields  and  rigid  body 
motions  are  evaluated  exactly,  which  accounts  for  the  terms  in  addition  to 
h  in  the  x  projection  vector,  see  Appendix  A.  However,  very  good  convergence 
can  be  obtained  in  many  cases  when  the  additional  terms  are  omitted.  Table  3 
compares  the  results  for  a  consistent  and  inconsistent  method  of  spurious  mode 
control  in  the  hemispherical  shell  problem;  the  inconsistent  method,  in  fact, 
converges  faster.  The  major  flaw  we  have  detected  in  inconsistent  methods  of 
stabilization  is  that  it  can  lead  to  erroneous  stresses  in  unusually  shaped 
elements  where  the  terms  anq  are  large. 
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Remark  6.2 : .  The  9-node  element  solutions  which  are  reported  here  for  the 
pinched  cylinder  with  diaphragm  could  be  obtained  by  simply  using  2x2 
quadrature.  When  the  assembled  stiffness  matrix  obtained  with  2x2 
quadrature  is  not  singular,  the  y-method  results  are  almost  identical  to  the 
results  obtained  by  2x2  quadrature. 


APPENDIX  A 


The  main  purpose  of  this  Appendix  is  to  provide  details  of  the  9-node 
shell  element  with  2x2  quadrature  and  the  stabilization  of  the  spurfous  modes 
by  the  y -method.  Stabilization  procedures  of  this  type  were  developed  for  the 
4-node  plate  by  8elytschko  and  Tsay  (1983).  Results  with  this 
y -element  have  been  reported  in  Belytschko  et  al .  (1985a)  but  the  consistency 
of  this  operator  for  curved  shells  was  not  examined  there;  formal  consistency 
is  equivalent  to  satisfying  the  patch  test,  see  Belytschko,  et  al .  (1984a). 

It  will  be  seen  here  that  once  the  ^-method  is  applied  to  curved  elements, 
certain  losses  of  consistency  result,  and  the  stabilization  can  best  be  called 
quasi -consi stent .  The  effect  of  this  is  not  clear,  but  as  can  be  seen  from 
the  results  presented  in  Section  6  no  deleterious  effects  have  been  noted. 

This  evidence  is  however  not  conclusive  since  even  inconsistent  control  can 
perform  well,  see  Table  4. 

The  degenerated  shell  element  is  developed  in  a  manner  that  closely 
follows  Hughes  and  Liu  (1981).  Attention  is  restricted  to  the  9-node 
element.  The  physical  domain  of  the  shell  is  described  by  an  isoparametric 
transformation  from  a  reference  cube  described  by  coordinates  £,  n,  C  through 
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X-j  =  I  h  ,  c)  ^-j  +  4*  j  ( C  •  n  ,  t)  (A.l) 


where  xi  and  x^_  are  the  coordinates  of  the  continuum  nodes  associated  with 

shell  node  I  and  x.  are  the  Cartesian  coordinates  in  the  physical  space;  see 


Fig.  Al .  As  can  be  seen  from  Fig.  Al ,  the  nodes  labeled  I-  are  below  the 
corresponding  shell  node  I  (;  <  0)  while  those  labeled  1+  above  the  shell  node 


I  (c  >  0). 
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The  shape  functions  ^ j  are  constructed  from  the  two-dimensional  Lagrange 
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interpolant  by 


^1+  s  N+(c)  N  j ( C  *  n) 
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ij(5)  -  j  5(5  -  1) 


<£,(0  » i  -  c2 
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(A. 2a) 


(A. 2b) 


(A. 2c) 


(A. 2d) 


(A.2e) 


(A. 3a) 


(A. 3b) 


(A. 3c) 


The  midsurface  of  the  shell  is  given  by  (A.l)  with  c  3  0,  which  by  means  of 
(A. 1-2)  can  be  written 
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xi I  =  2  (xil-  +  xi 1+^  ( A ,4b ) 

The  stiffness  matrix  will  be  integrated  numerically  and  the  points  at 
which  the  integrand  is  to  be  evaluated  are  called  quadrature  points.  At  each 

AAA 

quadrature  point,  a  local  Cartesian  coordinate  system  (x,  y,  z)  is  constructed 
so  that  the  x,  y  plane  is  tangent  to  the  midsurface  defined  by  (A. 4a)  and  the 
z  coordinate  is  perpendicular  to  this  plane.  If  large  rotations  of  the  shell 
are  to  be  treated,  this  coordinate  system  can  be  considered  corotational  in 
the  sense  used  by  Belytschko  and  Hsieh  (1973)  and  Belytschko  and  Marchertas 
(19/4),  but  this  is  not  of  major  concern  here. 

The  strain-displacement  equations  at  each  quadrature  point  can  be  written 
as 
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The  strain  matrix  is  defined  by 


*  T  r  An  *  *  *  *>  ^  . 

£  [£x»  ey>  2exy*  KX’  Ky ’  2lcxy’  exz’  eyz^  ^A’9^ 

A 

and  the  B  matrix  at  a  quadrature  point  a  gives  the  strains  through 
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ST  =  fi  (x  )d  =  l  Bj(x  )  dj  (A. 10) 

a  1*1  1  a  1 

d  a  ^2 ’  *  *  *  *  ig]  (A. 11a) 

I  “  ^uxl’  uyl *  uzl’  ♦xl*  ^yl^  (A. 11b) 


It  should  be  stressed  that  the  superposed  circumflex  in  the  above  equations 
indicates  that  all  components  of  e  and  d  are  expressed  in  the  corotational 
coordinate  system  of  point  a  when  evaluating  the  strain  at  point  a. 
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For  implementation  purposes,  a  transformation  R  is  defined  at  each 
quadrature  point  so  that 

ii  a  Si  (a. 12) 

so  that 

£  ■  j  *  I  1“  Si  Si  Ca.is) 

where  JR®  is  the  transformation  between  the  local  coordinate  system  of 
quadrature  point  a  and  the  coordinate  system  of  node  I.  If  5  d.o.f.  are  used 
per  node,  the  transformation  between  a  and  each  node  of  the  shell  will  be 
different  for  a  curved  element;  5  d.o.f.  per  node  is  necessary  if  singular 
systems  are  to  be  avoided. 

A 

The  B  matrix  is  given  by  the  component  matrices 


In  uniform  reduced  quadrature,  the  element  stiffness  is  given  by 


S,  ■  I  (E“)TlT<ia)  S(s.)  S(s„)  8°  •)<*„)  (A-ia) 

a«l 

I  , 

where  0  is  the  material  matrix  and  J  the  Jacobian  of  the  mapping  between  the 
reference  and  physical  volumes. 

Control  of  Spurious  Modes  for  Plate  Bending 

The  control  of  spurious  modes  for  the  9  node  plate  C°  element  has  been 
described  by  Belytschko,  et  al .  (1984b).  However,  the  method  described  there 
is  not  readily  applicable  to  the  curved  shell  element. 

In  order  to  describe  the  spurious  modes  and  their  control,  it  is 
convenient  to  define  row  vectors  b^,  the  elements  of  which  are  given  by 

!  blal  H  bxal  3  3V*a)/3*  (A*19a) 
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b2al  *  byal  =  3V~a)/3y 


(A. 19b) 


and  the  vectors  s  and  h  given  by 

sT  *  [1,  1,  1,  1,  1,  1,  1,  1,  1]  (A. 20) 

hT  *  [+1,  -1,  +1,  -1,  0,  -1,  +1,  -1,  +1]  (A. 21 ) 

A  single  coordinate  system  is  used  at  all  quadrature  points  of  the  plate 
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element  so  x  is  coincident  with  x,  y  with  y.  It  can  be  shown  that  if  the 
quadrature  points  are  xa,  which  in  the  reference  plane  are  given  by  the  four 
combinations  of  the  coordinates,  £  =>  n  3  ±  3'  ^  then 


bl  X  .  a  5  .  . 

~ia  ~J  ij 

(A. 22) 

bj  h  =  0 

~*1a  ~ 

(A. 23) 

bl  s  =  0 
-la  ~ 

(A. 24) 

N  h  *  N(x  )  h  »  -  i 

~<x  —  —  '-a  j 

(A. 25) 

N  S  »  N(x  )  S  *  1 
-a  ~  *** 

(A. 26) 

where  N  is  the  row  matrix  of  shape  functions,  5^  .  is  the  Kronecker  delta, 
x.  are  the  vectors  of  nodal  coordinates  of  the  element.  Using  (A. 22-26),  it 
can  easily  be  verified  that  the  first  3  modes  listed  in  Table  A1  are  spurious 
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singular  modes  of  the  plate  element,  which  means  that 


8  (x  )  d  ■  0 


(A. 27a) 


B  (x  )  d  =  0 


(A. 27b) 
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for  these  nodal  displacements  which  are  not  rigid  body  modes.  Note  that 
R“  *  ,1  for  a  plate,  and  that  {A. 25-26)  are  essential  in  establishing  (A. 27b). 

The  control  of  the  spurious  modes  then  consists  of  defining  a  projection 
operator  i  which  does  not  destroy  the  original  consistency  of  the  B  matrix,  as 
exhibited  by  (A. 22).  Following  Belytschko,  et  al .  (1984),  the  vector  j  is 

Q 

expanded  in  a  complete  representation  of  the  vector  space  R  through 


Table  Al 


Six  Spurious  Modes  of  the  9-Node  Element 
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(A. 28) 


Ten  vectors  are  needed  because  b.j  are  not  all  linearly  independent.  Linear 
consistency  then  requires  that 


XT(C0  s  +  c1  x  +  c2  x)  -  0 


for  all  c. 


(A. 29a) 


Substituting  (A. 28)  into  (A. 29a)  and  using  (A. 22-24)  yields 


co(al0  sT  s)  ♦  tj  (a,  hT  I  *la) 

a=l 


T  T  ^ 

c2  (a,  i  X  *  *io  i  i*  Z.  *2J  '  0 


(A. 29b) 


Since  the  above  must  hold  for  all  c.  ,  it  follows  that  a^Q  »  0  and  only  2 
conditions  need  be  satisfied  by  a^.  Thus  a  large  variety  of  options  are 
available  in  the  consistent  control  of  the  spurious  modes.  In  Belytschko,  et 
al .  (1984),  it  was  chosen  to  let 


ala  ■  -  T  *9  iT  s  a2a  "  •  T  a9  ST  *■  a  ■  1  t0  4  (A-30) 


However,  this  approach  is  not  adaptable  to  the  shell  problem,  as  will  become 
clear  shortly. 

The  approach  that  is  therefore  taken  is  to  define  a  projection 
operator  at  each  quadrature  point.  Equation  (A. 29)  must  hold  for  each 
and  this  is  accomplished  by  letting 
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3ia  = 


if  a  *  8 


-  ag  h  xi  if  a  =  S 


(A. 31) 


Hence 


x6  3  k  -  (iiT  £)  5xS  -  Qit  X)  ^ 


(A. 32) 


The  common  constant  a<j  has  been  omitted. 

In  this  scheme,  12  generalized  strains  (in  contrast  to  3  in  the  original 
scheme)  are  then  defined  by 
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where  3  identifies  the  quadrature  point. 
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Spurious  Mode  Control  for  9  Node  Shell  The  six  modes  which  are  to-  be 
controlled  are  listed  in  Table  Al.  The  form  of  a  spurious  mode  projection  has 
been  presented  in  Belytschko,  et  al .  (1985)  but  a  careful  examination  of  its 
consistency  and  limitations  is  still  not  available.  We  will  attempt  to 
partially  remedy  these  shortcomings  here. 

The  issues  of  consistency  (and  the  related  issues  of  patch  tests)  for 
shells  still  appear  unresolved.  While  the  requirement  that  an  element  be 
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strain-free  in  rigid  body  motion  is  quite  clear,  the  additional  consistency 
requirements  associated  with  linear  fields  are  more  difficult  to  formulate  in 
a  curved  element. 


For  this  reason,  we  have  taken  advantage  of  the  relationship  between  the 
shell  and  continuum  elements  (see  Fig.  Al)  to  develop  as  consistent  an 
operator  as  possible.  As  in  the  new  plate  version  of  the  spurious  mode 
projection  just  described,  the  projection  operator  is  defined  at  each 
quadrature  point  a,  and  denoted  by  ^  .  Consistency  then  requires  that 

T  3 

Ia  (co  2  +  l  cj  Xj)  •  0  for  all  c.  (A. 34) 
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where 


a„h  +  a  b  +  a  b 
9  -  xa  ~xa  ya  ~ya 


(A. 35) 


Note  that  we  have  immediately  dropped  the  big  matrices  associated  with 
quadrature  points  other  than  a,  as  in  Eq.  (A. 31),  and  a^,  which  must 
obviously  vanish  as  before.  Furthermore,  in  contrast  to  (A. 29),  consistency 
is  now  required  with  respect  to  fields  linear  in  z  (that  is  x^). 

We  now  define 
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(A. 36) 


that  is,  g1a  are  the  derivatives  of  the  continuum  node  shape  functions.  By 
consistency  of  the  original  continuum  shape  functions 


£ia  *j  =  >-SiaI+  Xjl  +  +  Sial-  Xjl-J  =  5i 


i,  j  *  1  to  3 


(A. 37) 


where  the  sum  is  over  the  continuum  nodes  and  r.  =  fx.  ,  x..),  which  are  the 

"J  ~j-  ~J  + 

coordinates  of  the  continuum  nodes. 

Using  (A. 2),  it  can  be  seen  that 
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for  i  =  1,  2  (A. 38) 


A  similar  expression  applies  to  the  nodes  below  the  midplane,  I-.  It  has  been 

A  A  A 

assumed  that  c  is  normal  to  the  (x,  y)  plane  so  that  3N+/3xi  =  0.  Since  this 
cannot  hold  exactly  the  approximate  equality  is  used  in  (A. 38). 

Using  (A. 38)  and  (A. 4b),  it  follows  that  for  i  *  1,  2;  j  *  1  to  3 
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Hence  using  (A. 19),  (A. 37)  and  (A. 39) 


Sia  *  5ij 


for  i  *  1,  2  j  =  1  to  3 


(A. 40) 


Employing  (A. 34-35)  with  (A. 40),  it  follows  that  the  requirements  for 
consistency  are 


Equations  (A.41a-b)  can  be  satisfied  exactly  as  in  the  development  of  the  flat 
plate  operator,  (A. 32),  so  that  the  projection  operator  at  each  point  becomes 


ia  3  5  -  (£T  *)  ixa  -  l)  5yS  (A-42) 

However,  (A. 41c)  cannot  be  satisfied  and  there  is  no  latitude  with  which  to 
attack  it,  so  that  the  method  as  described  here  must  violate  consistency 
whenever  z  *  0. 

The  control  of  the  spurious  modes  is  then  accomplished  by  defining  5 
additional  generalized  strains  at  each  of  the  four  laminar  quadrature  points 
through 
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The  first  two  generalized  strains  control  the  spurious  membrane  modes,  the 
next  three  control  the  bending  modes. 


The  stiffness  matrix  of  the  element  is  then  given  by  adding  the  effects 


of  the  stabilization  to  Ke  see  (A. 18),  which  gives 
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where  J  is  the  Jacobian  and  ra  is  a  diagonal  matrix  consisting  of 
blocks  r“  given  by 
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(A. 45) 

It  has  been  shown  by  Belytschko,  et  al .  (1984)  via  the  Hu-Washizu 

principle  that  the  y-method  stabilization  is  equivalent  to  recovering  that 

portion  of  the  diagonal  terms  of  the  stiffness  matrix  which  has  been  lost  by 

reduced  quadrature.  Specific  equations  for  the  r-parameters  have  been  given 

by  Belytschko,  et  al .  (1985a).  One  important  observation  on  the  choice  of  the 

r-parameters  is  that  while  ra  and  r  can  be  chosen  so  that  the  fully 
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integrated  stiffness  is  recovered,  r^  must  be  scaled  down  to  about  0.1; 
otherwise,  the  element  locks. 
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APPENDIX  B 


AN  EQUIVALENT  MIXED  ELEMENT 

We  show  here  a  construction  of  the  strain  and  stress  fields  such  that  the 
resulting  9-node  mixed  element  is  formally  equivalent  to  the  9-node  2x2 
quadrature  displacement  element.  The  demonstration  is  limited  to  the  membrane 
portion  of  the  9-node  element,  but  it  can  easily  be  extended  to  the  complete 
shell  element.  The  demonstration  is  the  spirit  of  the  Mai kus-Hughes  (1978) 
equivalence  theorem  but  offers  the  interesting  feature  that  by  using  Dirac 
functions  for  the  stress  interpolants,  the  demonstration  provides  an 
equivalence  with  the  exactly  integrated  mixed  element.  We  assume  that  the 
material  properties  of  the  element  are  constant. 

The  essential  ingredient  in  the  demonstration  of  the  equivalence  is  the 
construction  of  the  interpolation  functions  for  the  strains  and  stresses  which 
are  used  in  the  Hu-Washizu  functional.  For  conceptual  simplicity,  it  is 
important  to  choose  the  functions  so  that  the  orthogonality  condition  (3.8)  is 
met.  This  is  accomplished  as  follows. 

The  element  is  subdivided  into  4  subdomains  to  so  that  flj  contains 
integration  point  x^  and  no  two  subdomains  overlap,  i.e.  C\  flj  =  0  if 
I  *  J;  see  Fig.  B1  for  the  subdivision. 

We  then  define  Heaviside  functions  by 


1  if  x  is  In 

0  if  x  is  not  inn. 


(B.l) 
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The  Dirac  delta  function  is  denoted  by  <5(x  -  x^J.  Note  that 


/  Hj(x)  5(x  -  Xj)dn  =  I XJ  (B.2) 


where  IT1  are  the  components  of  the  unit  matrix  I. 

i  j  ~ 

The  strains  and  stresses  are  now.  interpolated  as  follows  (see  Fig.  81  for 
an  illustration  of  the  shape  functions) 


*  [£]_»  £2*  ~2* 


=  [S1#  Sjt  S3.  s4. 


where 


h  3  ¥*>  ~3x3  S  0  HI 


(B.3) 
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It  is  easily  verified  from  (3.8a)  and  (8,2)  that 


ri  =  /  ST  E.  dn  =  I . .  I , _ 

IJ  ~I  IJ  ~(3x3) 


or  5  *  1(12x12) 


From  (3.8b)  and  (8.6)  it  follows  that 


:^S 


(B.8) 


V  V  , 


VAV 
P" 


7^3 


•,v,\ 
■>  *. 


v^viN. 


>  .•-  .v\  .v  . 

.'■'.V.vV.V^ 


(B.9) 


m... 


r.  »  .  r.  v\ 


B  =  /  ST  B  da 
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from  (3.6c)  and  (B.5)  that 

‘a.d  0  0 

A.D  0 

_  2~  ~ 

~  symmetric  A^D 


where  Aj  are  the  areas  of  the  subdomains  a.. 

From  (3.12),  (B.9)  and  (B.10)  it  follows  that 


( B .10) 


K  *  BT  0  B 
/•»/<* 


(B.ll) 


which  is  analogous  to  the  form  of  the  numerically  integrated  displacement 
stiffness  except  that  Aj  replaces  the  determinant  of  the  Jacobian 
JI  3  J(i[b  see  (A. 18).  Note  that  the  weights  for  2x2  Gauss  quadrature  are 
all  unity.  However,  the  completion  of  the  proof  requires  that  it  be  shown 
that  flj  can  be  constructed  so  that  Jj  *  A^. 

To  show  that  for  an  arbitrary  isoparametric  quadratic  Lagrange  element 
the  partition  of  Q  into  flj  such  that  Aj  =  J^  is  possible,  we  observe  that 


J(C,  n)  =  det 


x. 


e 


{ B . 12a) 


is  a  cubic  function  of  both  5  and  n.  Thus  2x2  quadrature  evaluates  the  area 
of  the  element  exactly 

1  1  4 

A  3  /  /  J(e,  n)  dc  dn  3  [  J,  (B.12b) 

-1  -1  1=1  1 

The  above  equation  indicates  that  letting  A.  3  J.  we  have  l  A.  3  A.  The 

1  I  1 

partition  of  into  four  such  nonoverlapping  Qj  should  be  possible,  since  the 
only  restriction  is  that  each  integration  point  Xj  falls  within  sip 

It  is  also  of  interest  to  write  the  intermediate  equations,  namely,  the 
counterparts  of  (3.10).  Using  (B.7-10)  these  can  be  written  as 

er  3  B(xr)  d  (B.13a) 

Sj  3  Aj  D  ez  (B.13b) 

nE 

X  3  BT  s  3  l  BT(xr)  s.  (B.13c) 

1=1  1 

The  above  clearly  brings  out  the  structure  of  this  mixed  method.  The  strain 
parameters  are  simply  the  strains  at  the  quadrature  points,  the  stresses  are 
only  modified  by  the  areas  (Jacobi ans). 

An  almost  equivalent  mixed  method  can  be  constructed  by  letting 

h  *  ¥*^3x3  (B*14> 

Then  E  is  block-diagonal  with  the  blocks  given  by 


-II  "I  -(3x3) 
and  (B.9)  is  then  replaced  by 

St  *!> 

S<*2) 

S<s3> 

S(i4) 

where  by  the  mean  value  theorem,  is  some  point  in  ft.  (usually  not  the  Gauss 
quadrature  point).  Equations  analogous  to  (B.ll)  and  (B.13)  can  be  developed 
for  this  construction. 

It  is  not  clear  whether  this  construction  falls  within  the  scope  of  the 
convergence  proof  given  by  Oden  and  Reddy  (1974).  However,  the  ucs  of  Dirac 
delta  functions  makes  this  method  a  hybrid  of  collocation  and  Galerkin 
methods,  for  which  convergence  theory  should  be  feasible. 
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APPENDIX  C 

C.l  Impulsively  Loaded  Cylindrical  Panel 

This  problem  is  defined  in  Fig.  C.l  and  Table  C.2.  The  Ends  of  the  panel 
are  simply-supported,  while  the  side  boundaries  are  fixed.  Both  material  and 
geometric  nonlinearities  were  considered.  The  load  was  applied  by  prescribing 
the  initial  velocity  given  in  Fig.  C.2  to  the  nodes  in  the  region  loaded  by 
the  explosive.  Figure  C.3  shows  an  undeformed  and  deformed  mesh.  Table  C.2 
compares  the  results  obtained  for  various  meshes  to  an  experimental  result. 

It  can  be  seen  that  convergence  to  the  experimental  value  is  relatively 
slow.  The  reason  for  this  however,  is  not  the  accuracy  of  the  element,  but 
tne  extreme  localization  of  deformation  which  occurs  due  to  the  formation  of 
plastic  hinges. 

C.2  Collapse  of  a  Hollow  Column 

Figure  C.4  shows  the  simulation  of  a  hollow  column  loaded  axially.  The 
time  history  of  the  load  and  material  and  geometric  properties  may  be  found  in 
Fig.  3  and  Table  1  of  Kennedy  and  Belytschko  (1982)  where  results  obtained 
with  a  4-node  quadrilateral  element  using  one-point  quadrature  and  hourglass 
control  are  also  given.  The  results  obtained  with  this  element  compare  well 
with  those  obtained  for  the  quadrilateral,  except  the  model  is  somewhat 
stiffer.  Note  the  severe  change  n  cross-section  which  accompanies  buckling. 
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TABLE  C.1 

Material  and  Geometric  Constants  for  Impulsively 
Loaded  Cylindrical  Panel 


Young's  modulus 
Density 

Poisson's  ratio 
Yield  stress 
Plastic  modulus 
Impulse  over  Rj 


E  =  10.5  x  106 

p  =  2.5  x  104  lb-sec2/in4 

v  *  0.33 

a  *  44000.  psi 

E  *  0.  psi 
p 

5650  in/sec 


TABLE  C. 2 


Final  Displacements  of  the  Cylindrical  Panel 
for  Various  Meshes  with  Ilyushin  Yield  Condition 


M 

Mesh 

Displacement 

Displacement 

CPU 

time 

P 

Hal f-panel 

at  y  =  6.28 

at  y  -  9.42 

IBM 

3033 

1 

6  x  16 

0.917 

0.401 

83 

Sa 

8  x  16 

1.043 

0.448 

138 

£ 

10  x  20 

1.081 

0.462 

162 

h/. 

12  x  24 

1.124 

0.473 

291 

'r* 

16  x  32 

experimental  [48] 

1.183 

1.28 

0.530 

0.70 

670 
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